{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b0 = -67.3538\n",
      "b12=106.9354    0.0275\n",
      "拟合优度R2=0.9836\n"
     ]
    }
   ],
   "source": [
    "#12.1\n",
    "#使用线性回归解决\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression\n",
    "a = np.array([\n",
    "    [4,3424,450.5],\n",
    "    [5,4086,613.9],\n",
    "    [4,4388,501.5],\n",
    "    [7,4808,781.5],\n",
    "    [5,5896,611.1],\n",
    "    [10,6604,1222.1],\n",
    "    [7,6662,793.2],\n",
    "    [6,7018,792.7],\n",
    "    [9,8706,1121],\n",
    "    [8,10478,1094.2]])\n",
    "model = LinearRegression().fit(a[:,:2],a[:,2])\n",
    "y = model.predict(np.array([10,9600]).reshape(1,-1))\n",
    "b0 = model.intercept_ #常数系数\n",
    "b12 = model.coef_#k系数\n",
    "R2 = model.score(a[:,:2],a[:,2])#计算R2\n",
    "print('b0 = %.4f\\nb12=%.4f%10.4f'%(b0,b12[0],b12[1]))\n",
    "print('拟合优度R2=%.4f'%R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "b0 = 62650.9285\n",
      "b12=-979.0570    0.2862\n",
      "拟合优度R2=0.9022\n"
     ]
    }
   ],
   "source": [
    "#12.1\n",
    "\n",
    "import numpy as np\n",
    "from sklearn.linear_model import LinearRegression\n",
    "a = np.array([\n",
    "    [23.56,76200,59190],\n",
    "    [24.44,91200,65450],\n",
    "    [32.07,106700,62360],\n",
    "    [32.46,111600,64700],\n",
    "    [31.15,119000,67400],\n",
    "    [34.14,129200,64440],\n",
    "    [35.30,143400,68000],\n",
    "    [38.7,159600,72400],\n",
    "    [39.63,180000,75710],\n",
    "    [46.68,193000,70680]])\n",
    "model = LinearRegression().fit(a[:,:2],a[:,2])\n",
    "y = model.predict(np.array([10,9600]).reshape(1,-1))\n",
    "b0 = model.intercept_ #常数系数\n",
    "b12 = model.coef_#k系数\n",
    "R2 = model.score(a[:,:2],a[:,2])#计算R2\n",
    "print('b0 = %.4f\\nb12=%.4f%10.4f'%(b0,b12[0],b12[1]))\n",
    "print('拟合优度R2=%.4f'%R2)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 432x288 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEGCAYAAABmXi5tAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAATI0lEQVR4nO3df4xVaX3H8c+Hu0xpB1DTGaebnVmZmCGF+gNxoNts0mykm8VtA02cAEttu02VtILaamrXxqyITZP2D2kqFDNVW2txcbM1LiqGJXFbUtEts4goP1Z+aZeFzc6sulQ6usOdb/+4l/XO7B24wz33njvzvF/JJPc858l9vg8HPnM4Px0RAgDMfnPyLgAA0BwEPgAkgsAHgEQQ+ACQCAIfABJxS94FTKWjoyMWLVqUdxkAMKM8+eSTIxHRWW1dywb+okWLNDQ0lHcZADCj2P7BVOs4pAMAiSDwASARBD4AJILAB4BEEPgAkAgCH0AuiuOhwYNntWzbYxo8eFbFcR7k2Ggte1kmgNnr/MgVbd59ROdHrmh0rKjtB07r0aMXtWPjcvV2tOdd3qzFHj6AphvYdUinnr2s0bGiJGl0rKiTly5rYNehnCub3Qh8AE3X1zVfk4/gjIe0uGtBPgUlIpPAt73a9lO2z9h+oMr6220/bvtbto/ZvjeLcQHMTOtX9Ki9rTChrb2toHUrunOqKA11B77tgqSdkt4qaamk+2wvndTtQ5Iejog3Sdog6R/rHRfAzLVqSZcKczyhrTDHWrWkK6eK0pDFSduVks5ExDlJsr1H0lpJJyr6hKSF5c+vkHQxg3EBzFAL583Vsa335F1GcrI4pHObpKcrli+U2yptlfR22xck7ZP07mpfZHuT7SHbQ8PDwxmUBgC4plknbe+T9C8R0S3pXkmftf2ysSNiMCL6I6K/s7Pq0z0BYNZq9L0JWRzSeUZST8Vyd7mt0h9LWi1JEfEN2/MkdUh6LoPxAWDGa8a9CVns4R+W1Ge713abSidl907q8z+SVkmS7SWS5knimA0AlDXj3oS6Az8irkraImm/pJMqXY1z3PY222vK3d4v6Z22vy3pIUn3RwT3UQNAWTPuTcjk0QoRsU+lk7GVbQ9WfD4h6c4sxgKA2Wj9ih5958ILuvJi8aW2rO9N4E5bAGgBzbg3gYenAUALaMa9CezhA0AiCHwASASBDwCJIPABIBEEPgAkgsAHgEQQ+ACQCAIfABJB4ANAIgh8AEgEgQ8AiSDwASARBD4AJILAB4BEEPgAkAgCHwASQeADQCIIfABIBIEPAIkg8AEgEQQ+ACSCwAeARBD4AJAIAh8AEkHgA0AiCHwASASBDwCJIPABIBEEPgAkIpPAt73a9lO2z9h+oMr67baPln++Z/vHWYwLAKjdLfV+ge2CpJ2S7pZ0QdJh23sj4sS1PhHx5xX93y3pTfWOCwCYniz28FdKOhMR5yLiRUl7JK29Tv/7JD2UwbgAgGnIIvBvk/R0xfKFctvL2H6NpF5JX5ti/SbbQ7aHhoeHMygNAHBNs0/abpD0SEQUq62MiMGI6I+I/s7OziaXBgCzWxaB/4yknorl7nJbNRvE4RwAyEUWgX9YUp/tXtttKoX63smdbP+qpFdJ+kYGYwIApqnuwI+Iq5K2SNov6aSkhyPiuO1tttdUdN0gaU9ERL1jAgCmr+7LMiUpIvZJ2jep7cFJy1uzGAsAcHO40xYAEkHgA0AiCHwASASBDySuOB4aPHhWy7Y9psGDZ1Uc57qK2SqTk7YAZqbzI1e0efcRnR+5otGxorYfOK1Hj17Ujo3L1dvRnnd5yBh7+EDCBnYd0qlnL2t0rHTz++hYUScvXdbArkM5V4ZGIPCBhPV1zdfkIzjjIS3uWpBPQWgoAh9I2PoVPWpvK0xoa28raN2K7pwqQiMR+EDCVi3pUmGOJ7QV5lirlnTlVBEaiZO2QMIWzpurY1vvybsMNAl7+ACQCAIfABJB4ANAIgh8AEgEgQ8AiSDwASARBD4AJILAB4BEEPgAkAgCHwASQeADQCIIfABIBIEPAIkg8AEgEQQ+ACSCwAeARBD4AJAIAh8AEkHgA0AiCHwASEQmgW97te2nbJ+x/cAUfdbZPmH7uO3PZTEucDOK46HBg2e1bNtjGjx4VsXxyLskoCluqfcLbBck7ZR0t6QLkg7b3hsRJyr69En6oKQ7I+JHtl9d77jAzTg/ckWbdx/R+ZErGh0ravuB03r06EXt2LhcvR3teZcHNFQWe/grJZ2JiHMR8aKkPZLWTurzTkk7I+JHkhQRz2UwLjBtA7sO6dSzlzU6VpQkjY4VdfLSZQ3sOpRzZUDjZRH4t0l6umL5Qrmt0mJJi21/3fY3ba+u9kW2N9kesj00PDycQWnARH1d8zX5CM54SIu7FuRTENBEzTppe4ukPkl3SbpP0j/ZfuXkThExGBH9EdHf2dnZpNKQkvUretTeVpjQ1t5W0LoV3TlVBDRPFoH/jKSeiuXuclulC5L2RsRYRJyX9D2VfgEATbVqSZcKczyhrTDHWrWkK6eKgOap+6StpMOS+mz3qhT0GyRtnNTniyrt2f+z7Q6VDvGcy2BsYFoWzpurY1vvybsMIBd17+FHxFVJWyTtl3RS0sMRcdz2Nttryt32S3re9glJj0v6i4h4vt6xAQC1c0RrXoPc398fQ0NDeZcBADOK7Scjor/aOu60BYBE1HQM3/ZDKv1y6JXUJeldEfGVRhYGAMhWrXv4b5R0LiJWSvo9SR9uXEkAgEa4YeDbniepU9JHyk0nJL2qkUUBALJXyx7+6ySdjoiflpeXS/p240oCADRCLcfw3yjp9vKefkGlPf0PNLQqAEDmag38L0h6QtJcSX8TEV9vaFUAgMzVGvibIuI9jS4GANA4tRzDf62k040uBADQWDfcw48IHiMIALMAd9oCQCIIfABIBIEPAIkg8AEgEQQ+ACSCwAeARBD4AJAIAh8AEkHgA0AiCHwASASBDwCJIPABIBEEPgAkgsAHgEQQ+ACQCAIfABJB4ANAIgh8AEgEgQ8AiSDwASARBD4AJCKTwLe92vZTts/YfqDK+vttD9s+Wv55RxbjAgBqd0u9X2C7IGmnpLslXZB02PbeiDgxqevnI2JLveMBAG5OFnv4KyWdiYhzEfGipD2S1mbwvQCADGUR+LdJerpi+UK5bbK32T5m+xHbPdW+yPYm20O2h4aHhzMoDQBwTbNO2n5J0qKIeIOkA5I+U61TRAxGRH9E9Hd2djapNABIQxaB/4ykyj327nLbSyLi+Yj4WXnxk5LenMG4AIBpyCLwD0vqs91ru03SBkl7KzvYvrVicY2kkxmMCwCYhrqv0omIq7a3SNovqSDp0xFx3PY2SUMRsVfSe2yvkXRV0g8l3V/vuACA6XFE5F1DVf39/TE0NJR3GQAwo9h+MiL6q63jTlsASASBDwCJIPABIBEEPgAkgsAHgEQQ+ACQCAIfABJB4ANAIgh8AEgEgQ8AiSDwASARBD4AJILAB4BEEPgAkAgCHwASQeBjWorjocGDZ7Vs22MaPHhWxfHWfJ8CgJer+41XSMf5kSvavPuIzo9c0ehYUdsPnNajRy9qx8bl6u1oz7s8ADfAHj5qNrDrkE49e1mjY0VJ0uhYUScvXdbArkM5VwagFgQ+atbXNV+Tj+CMh7S4a0E+BQGYFgIfNVu/okftbYUJbe1tBa1b0Z1TRQCmg8BHzVYt6VJhjie0FeZYq5Z05VQRgOngpC1qtnDeXB3bek/eZQC4SezhA0AiCHwASASBDwCJIPABIBEEPgAkgsAHgEQQ+ACQCAIfABKRSeDbXm37KdtnbD9wnX5vsx22+7MYFwBQu7oD33ZB0k5Jb5W0VNJ9tpdW6bdA0nslPVHvmACA6ctiD3+lpDMRcS4iXpS0R9LaKv0+KulvJf00gzEBANOUReDfJunpiuUL5baX2F4uqScivnK9L7K9yfaQ7aHh4eEMSgMAXNPwk7a250j6mKT336hvRAxGRH9E9Hd2dja6NABIShaB/4yknorl7nLbNQskvU7Sf9j+vqQ7JO3lxC0ANFcWgX9YUp/tXtttkjZI2nttZUS8EBEdEbEoIhZJ+qakNRExlMHYAIAa1R34EXFV0hZJ+yWdlPRwRBy3vc32mnq/HwCQjUxegBIR+yTtm9T24BR978piTADA9HCnLQAkgsAHgEQQ+ACQiFkV+MXx0ODBs1q27TENHjyr4njkXRIAtIxMTtq2gvMjV7R59xGdH7mi0bGith84rUePXtSOjcvV29Ged3kAkLtZs4c/sOuQTj17WaNjRUnS6FhRJy9d1sCuQzlXBgCtYdYEfl/XfE0+gjMe0uKuBfkUBAAtZtYE/voVPWpvK0xoa28raN2K7pwqAoDWMmsCf9WSLhXmeEJbYY61aklXThUBQGuZNSdtF86bq2Nb78m7DABoWbNmDx8AcH0EPgAkgsAHgEQQ+ACQCAIfABJB4ANAIgh8AEgEgQ8AiSDwASARBD4AJILAB4BEEPgAkAgCHwASQeADQCIIfABIBIEPAIkg8AEgEQT+dRTHQ4MHz2rZtsc0ePCsipPfkg4AM8isecVh1s6PXNHm3Ud0fuSKRseK2n7gtB49elE7Ni5Xb0d73uUBwLSxhz+FgV2HdOrZyxodK0qSRseKOnnpsgZ2Hcq5MgC4OQT+FPq65mvyEZzxkBZ3LcinIACoUyaBb3u17adsn7H9QJX1f2L7O7aP2v4v20uzGLeR1q/oUXtbYUJbe1tB61Z051QRANSn7sC3XZC0U9JbJS2VdF+VQP9cRLw+IpZJ+jtJH6t33EZbtaRLhTme0FaYY61a0pVTRQBQnyxO2q6UdCYizkmS7T2S1ko6ca1DRFyu6N8uqeUvd1k4b66Obb0n7zIAIDNZBP5tkp6uWL4g6dcnd7K9WdL7JLVJeku1L7K9SdImSbr99tszKA0AcE3TTtpGxM6IeK2kv5T0oSn6DEZEf0T0d3Z2Nqs0AEhCFoH/jKSeiuXucttU9kj63QzGBQBMQxaBf1hSn+1e222SNkjaW9nBdl/F4m9LOp3BuACAaaj7GH5EXLW9RdJ+SQVJn46I47a3SRqKiL2Sttj+LUljkn4k6Q/rHRcAMD2ZPFohIvZJ2jep7cGKz+/NYhwAwM1zRGteIWl7WNIPGjxMh6SRBo/RaMyhNTCH1jEb5lHPHF4TEVWvemnZwG8G20MR0Z93HfVgDq2BObSO2TCPRs2BZ+kAQCIIfABIROqBP5h3ARlgDq2BObSO2TCPhswh6WP4AJCS1PfwASAZBD4AJGLWB77tT9t+zvZ3p1hv2/9QfnnLMdvLm13jjdQwh7tsv1B+wcxR2w9W65cn2z22H7d9wvZx2y+7Ga/Vt0WNc2jpbWF7nu3/tv3t8hw+UqXPL9j+fHk7PGF7UQ6lTqnGOdxve7hiO7wjj1pvxHbB9rdsf7nKuuy3Q0TM6h9JvylpuaTvTrH+XklflWRJd0h6Iu+ab2IOd0n6ct513mAOt0paXv68QNL3JC2dSduixjm09LYo/9nOL3+eK+kJSXdM6vMuSZ8of94g6fN5130Tc7hf0o68a61hLu+T9Llqf2casR1m/R5+RByU9MPrdFkr6V+j5JuSXmn71uZUV5sa5tDyIuJSRBwpf/5fSSdVepdCpZbeFjXOoaWV/2x/Ul6cW/6ZfOXGWkmfKX9+RNIq21aLqHEOLc92t0oPk/zkFF0y3w6zPvBrUO0FLjPqH3HZb5T/i/tV27+WdzHXU/6v6ZtU2jOrNGO2xXXmILX4tigfRjgq6TlJByJiyu0QEVclvSDpl5ta5A3UMAdJelv50OAjtnuqrM/b30v6gKTxKdZnvh0I/NnhiErPz3ijpI9L+mK+5UzN9nxJ/y7pz2Liqy9njBvMoeW3RUQUo/R+6W5JK22/LueSpq2GOXxJ0qKIeIOkA/r5nnJLsP07kp6LiCebOS6BP/0XuLSciLh87b+4UXpy6VzbHTmX9TK256oUlLsj4gtVurT8trjRHGbKtpCkiPixpMclrZ606qXtYPsWSa+Q9HxTi6vRVHOIiOcj4mflxU9KenOTS7uROyWtsf19lV4K9Rbb/zapT+bbgcAvvazlD8pXiNwh6YWIuJR3UdNh+1euHduzvVKl7dpS/0DL9X1K0smI+NgU3Vp6W9Qyh1bfFrY7bb+y/PkXJd0t6dSkbnv183dWDEj6WpTPHLaCWuYw6dzPGpXOt7SMiPhgRHRHxCKVTsh+LSLePqlb5tshk+fhtzLbD6l05USH7QuSPqzSSR5FxCdUeo7/vZLOSPo/SX+UT6VTq2EOA5L+1PZVSaOSNrTSP9CyOyX9vqTvlI+9StJfSbpdmjHbopY5tPq2uFXSZ2wXVPpl9HBEfNkTX1j0KUmftX1GpYsFNuRXblW1zOE9ttdIuqrSHO7PrdppaPR24NEKAJAIDukAQCIIfABIBIEPAIkg8AEgEQQ+ACSCwAeARBD4AJAIAh+oUflZ+HeXP/+17Y/nXRMwHbP+TlsgQx+WtM32q1V6UuaanOsBpoU7bYFpsP2fkuZLuqv8THxgxuCQDlAj269X6TkuLxL2mIkIfKAG5acv7lbpLUQ/sT35kcJAyyPwgRuw/UuSviDp/RFxUtJHVTqeD8woHMMHgESwhw8AiSDwASARBD4AJILAB4BEEPgAkAgCHwASQeADQCL+Hz+3hw5imjWvAAAAAElFTkSuQmCC\n"
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   R-squared:                       0.929\n",
      "Model:                            OLS   Adj. R-squared:                  0.915\n",
      "Method:                 Least Squares   F-statistic:                     65.60\n",
      "Date:                Mon, 04 Jul 2022   Prob (F-statistic):           0.000465\n",
      "Time:                        23:27:31   Log-Likelihood:                 11.119\n",
      "No. Observations:                   7   AIC:                            -18.24\n",
      "Df Residuals:                       5   BIC:                            -18.35\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "const          0.1781      0.053      3.349      0.020       0.041       0.315\n",
      "x1             0.1673      0.021      8.100      0.000       0.114       0.220\n",
      "==============================================================================\n",
      "Omnibus:                          nan   Durbin-Watson:                   2.563\n",
      "Prob(Omnibus):                    nan   Jarque-Bera (JB):                0.683\n",
      "Skew:                           0.287   Prob(JB):                        0.711\n",
      "Kurtosis:                       1.581   Cond. No.                         6.99\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[-1.40066042  0.73900674] \n",
      " 剂量为： 1.8953283452155225\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\goodboy\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\statsmodels\\stats\\stattools.py:74: ValueWarning: omni_normtest is not valid with less than 8 observations; 7 samples were given.\n",
      "  warn(\"omni_normtest is not valid with less than 8 observations; %i \"\n"
     ]
    }
   ],
   "source": [
    "#12.3\n",
    "\n",
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "import pylab as plt\n",
    "a = np.array([[0.9,0.37],[1.1,0.31],[1.8,0.44],[2.3,0.60],[3.0,0.67],[3.3,0.81],[4.0,0.79]])\n",
    "x = a[:,0]; p = a[:,1]\n",
    "plt.plot(x, p, 'p');\n",
    "plt.xlabel('$x$')\n",
    "plt.ylabel('$p$', rotation=0);\n",
    "plt.show()\n",
    "X = sm.add_constant(x)\n",
    "md = sm.OLS(p, X).fit()  #构建并拟合模型\n",
    "print(md.summary())   #输出计算结果\n",
    "y = np.log(p/(1-p))\n",
    "cs = np.linalg.pinv(X) @ y  #线性最小二乘拟合\n",
    "x0 = -cs[0] / cs[1]   #求剂量\n",
    "print(cs,'\\n','剂量为：',x0)\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "拟合的参数为： [-1.09906293  0.28603066  1.78234146 -0.3619957 ]\n",
      "两种水分含量分别为： 8.994643947595016 , 3.842465455803507\n",
      "加盖和不加盖的赔率分别为： 1.255442704370899 , 1.853583938291554\n",
      "赔率比为： 0.6773055584027337\n"
     ]
    }
   ],
   "source": [
    "#12.4\n",
    "\n",
    "import numpy as np\n",
    "a = np.loadtxt('ti10_7.txt')      #提取10×8矩阵\n",
    "d = np.vstack([a[:,:4], a[:,4:]]) #变形为20×4矩阵\n",
    "f = d[:,3] * (d[:,2]==1) +(100-d[:,3])*(d[:,2]==0) #发芽频数\n",
    "pd = f/100  #发芽频率\n",
    "mat = np.stack([np.ones(20), d[:,0],d[:,1], d[:,0]*d[:,1]]).T\n",
    "cs = np.linalg.pinv(mat) @ np.log(pd/(1-pd))  #拟合参数\n",
    "print('拟合的参数为：', cs)\n",
    "x11 = -(cs[0]+cs[2]) / (cs[1]+cs[3])  #计算水分含量\n",
    "x10 = -cs[0] / cs[1]                  #计算水分含量\n",
    "print('两种水分含量分别为：', x11,',',x10)\n",
    "odd1 = np.exp(cs@[1,6,1,6]) #计算加盖时发芽的赔率\n",
    "odd2 = np.exp(cs@[1,6,0,0])  #计算不加盖时发芽的赔率\n",
    "OR = odd1/odd2              #计算赔率比\n",
    "print('加盖和不加盖的赔率分别为：',odd1,',',odd2)\n",
    "print('赔率比为：', OR)\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}